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ABSTRACT 

We explore the self-regulation of star formation using a large suite of high resolution hydrodynamic 
simulations, focusing on molecule-dominated regions (galactic centers and [U]LIRGS) where feedback 
from star formation drives highly supersonic turbulence. In equilibrium the total midplane pressure, 
dominated by turbulence, must balance the vertical weight of the ISM. Under self-regulation, the 
momentum flux injected by feedback evolves until it matches the vertical weight. We test this flux 
balance in simulations spanning a wide range of parameters, including surface density S, momentum 
injected per stellar mass formed (p*/m*), and angular velocity. The simulations are two dimensional 
radial- vertical slices, and include both self-gravity and an external potential that helps to confine gas 
to the disk midplane. After the simulations reach a steady state in all relevant quantities, including 
the star formation rate Esfr, there is remarkably good agreement between the vertical weight, the tur- 
bulent pressure, and the momentum injection rate from supernovae. Gas velocity dispersions and disk 
thicknesses increase with p*/m„. The efficiency of star formation per free-fall time at the mid-plane 
density, eg (no), is insensitive to the local conditions and to the star formation prescription in very 
dense gas. We measure eff(no)~0. 004-0. 01, consistent with low and approximately constant efficien- 
cies inferred from observations. For Eg(100-1000) M Q pc -2 , we find Esfr€(0.1-4) M kpc~ 2 yr _1 , 
generally following a Esfroc E 2 relationship. The measured relationships agree very well with verti- 
cal equilibrium and with turbulent energy replenishment by feedback within a vertical crossing time. 
These results, along with the observed E-Esfr relation in high density environments, provide strong 
evidence for the self-regulation of star formation. 

Subject headings: galaxies: ISM - galaxies: kinematics and dynamics - galaxies: starburst - galaxies: 
star formation - ISM: structure - turbulence 
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1. INTRODUCTION 

1.1. Star Formation on Galactic Scales 

Observations reveal that stars form in the molecular 
component of the interstellar medium (ISM). Therefore, 
the dynamics of molecular gas must affect the star forma- 
tion process. On galactic scales, gravity concentrates gas 
into clouds in which stars eventually form. The resulting 
feedback from stellar winds, ionizing and non-ionizing ra- 
diation, and supernovae (SN) (either local or nearby in 
the disk) redisperses this dense gas. The formation, de- 
struction, and the dynamical state of star forming clouds 
depend strongly on the local conditions of the ISM. In 
(ultra) luminous infrared galaxies ([U]LIRGs) and the 
centers of galaxies, molecular gas pervades much of the 
ISM, including regions not actively forming stars. Gas in 
such environments has higher mean volume and surface 
density compa red to the gas found in giant molecular 
clouds (GMCs, iSolomon fe Vanden Boutll2005f) in lower- 
density regions of galaxies. Near-future ALMA observa- 
tions will resolve high density tracers, and thereby reveal 
the detailed structure and kinematics of gas in starbursts. 
Understanding how small-scale feedback associated with 
star-formation acts in concert with larger scale processes 
in starbursts (as well as mid- and outer- disks) is cru- 
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cial for developing any successful theory of galactic star 
formation. 

Stellar feedback plays a key role in regulating the 
ther mal balance and morphological structure of th e 
ISM (|McKee fe Ostriker|[l977l: iNorman fc Lkeuchilll989ft . 
Feedback is also believe d to be the primary mec hanism 
driving turbulence (e.g. INorman fc F crrara 199(1). Since 
turbulence is observed on all scales larger than the size of 
the densest prestellar cores, it is now understood to be an 
essential component controlling the dynamics and regu- 
lating star formation in the ISM (see lMac Low fc Klessenl 
l20Q4TlMcKee fc Ostrikerl [20071 and references therein). 
The vertical scale height of the galactic disk depends 
on the balance between gaseous, stellar, and dark mat- 
ter potentials that concentrate gas, and the pressures 
(thermal, turbulent, magnetic, cosmic ray, and radia- 
tion) that oppose grav ity and limit runaway collapse (e.g. 
IBoulares fc C^ll990l) . 

Over sufficiently large scales, the star formation rate 
surface density, Esfr, is obse rved to correla t e well with 
the gas su rface density, E (|Schmidtl 119591 : iKennicutH 
119891 [19981) . This correlation appears to take on various 
forms in different regions within disk galaxies. In the 
outer-disk regions containing little molecular gas, there 
is no univers al power-law index describing the Esfr— E 
relationship (|Bigiel et al.ll2010l ). Ins tead, Esfr depends 
on both E and the stellar density (jBlitz fc Rosolowskvl 
12004 120061) : this is presumably because stellar rather 
than gas vertical gravity dominates in outer disks (see 
below). At smaller radii, by mass the ISM is domi- 
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nated by molecular gas, for which two different star for- 
mation laws appear to take hold. In mid-disk regions 
where most of the volume is filled with atomic gas and 
molecules are confine d in isolated GMCs ( with a limited 
range of properties - iSheth et al.1 (|2008l ): iBolatto et al.1 
(2008)), there is a strong, approximately linear correla- 
tion between th e star formation rate and molecular mass , 
Esfr. oc E mo j (IWong fc Blitd l200"i IBigiel et all [2008; 
ISchruba et al.ll201lD . Towards the central regions and 
in sta rbursts where the ISM is almos t completely molec- 
ular ([Solomon fc Vanden Boutl 120051 ) , there appears to 
be a steeper Esfr - E mo i relatio nship (IKennicuttl 119981: 
Genzel et all [2010: Da ddi et all 120101: iNaravanan et al.l 
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The variations in E-Esfr correlations in different 
galactic regions presumably owe to the differences in the 
characteristics of the ISM. Gas properties such as the 
temperatures, densities, and velocities are found to vary 
between starbursts and more quiescent environments. In 
the Gala ctic center, molecular g as is much more preva- 
lent (e.g. iBallv et al.|[l987l ll988!). and Esfr is measured 
to be ~1.5 dex higher t han in the mid- to outer- disk 
(|Yusef-Zadeh et al.H2009f) . Observed linewidths from the 
dense, molecular gas in the Ga lactic Center are mea- 
sured to reach £ 10 km s" 1 ([Oka et al.l [l99l [200H: 
IShettv et alll2012D. and as hi g h as ^ 100 km s" 1 in star- 
bursts (iSolomon et all 119971: iDownes fc Solomon! 119981: 
IGenzel et al.l 120111 ). Turbulent velocities in GMCs are 
signif i cantly lower, ranging from 1-6 km s -1 (e.g. lLarsonl 
119811 : [Solomon et al.l 119871) . However, present observa- 
tions of (U)LIRGs do not have sufficient resolution to 
distinguish between perturbed motions (such as large- 
scale streaming) on scales /t H, the disk thickness, and 
more localized turbulence (i.e. velocity dispersions on 
^10 pc scales, similar to GMCs). 

Global numerical simulations of disk galaxies have 
show n that the structures formed by self-gravity 
(e.g. IShettv fc Ostriker! 12001 iDobbs fc PringleT [20091 
ITasker fc Tanl 120091: iTaskerl 120111) or by cloud colli- 
sions (e.g. IDobbsl 120081: iTasker fc Tanl 120091 ) can gener- 
ally reproduce observed morphological features of the 
ISM, such as filamentary substructure, cloud masses, 
sizes, and basic kinematic properties. Additionally, 
large scale simulations have suggested that gravita- 
tional instability naturally results in power-law relation- 
ships between Esfr and E if the Toomre Q and ve- 
locity dispersion are un iform (e.g. ILi et al.l 120051 120061 : 
IShettv fc Ostrikerl 120081) . Simulations with feedback 
have produced a range in the exponent and coeffi- 
cient of the Esfr-E relationship, depending on the spe- 
cific f eedback prescription (e . g ITasker fc Brvanl 120061 
2008 1 : [Robertson fc Kravtsovl 120081: IShettv fc Ostrikerl 



20081 iKovama fc Ostrikerl I2009at [Pobbs et all 120111: 
Hopkins et all 1201 IShettv fc Ostrikerl (|2008l) pointed 



out that the relationship between Esfr an d E in gen- 
eral should depend on the thickness of the gas disk, and 
therefore on the gas velocity dispersion and on the stellar 
potential if it dominates (see below). 

Variations in feedback parameters, such as the injected 
momenta, energies, and rates, combined with other pro- 
cesses such as rotation, vertical motions due to an exter- 
nal potential, shear, and large-scale gravitational insta- 
bility in the shearing, rotating flow, are likely to con- 



tribute to the observed differences in velocity disper- 
sions between sta r bursts and more quie s cent regions. 
IQstriker fc Shettvl (|20ll and IKim et all (I20TI (here- 
after KKOll) argue that the velocity dispersion on scales 
comparable to the neutral gas disk's thickness will be 
relatively constant if turbulence is driven by feedback, 
because the driving rate and dissipation rate both scale 
inversely with the vertical crossing time (or gravitational 
free-fall time) of the ISM. Simulations of the ISM in mid- 
and outer- disk environments have shown that velocity 
dispersions are in fac t not strongly sensitive to the feed- 
back parameters (e.g. IDib et alj|2006t IShettv fc Ostrikerl 



120081: 1 Joung et all2009l KKOll). Such simulations allow 
for a detailed assessment of the relationships between the 
relevant physical quantities, and provide a direct avenue 
for testing analytical theories of star formation. 

1.2. Theory of Star Formation Self- Regulation 

A theory for the self-regulation of star forma- 
tion on galactic scal e s has recently been formu- 
lated bv IQstriker et all (120101) (hereafter OML10) and 
IQstriker fc Shettvl ([201 ll ) (hereafter Paper I). KKOll 
conducted numerical models of multi-phase gaseous disks 
in the regime where diffuse atomic gas dominates (X ,$ 
20 M Q pc -2 ), verifying the assumptions and predicted 
features of the self-regulated thermal/dynamical equilib- 
rium theory. In the present work, we shall instead fo- 
cus on numerical simulations of the molecule-dominated 
starburst regime. To provide an overall context and dis- 
tinguish between the various regimes, we briefly review 
the concepts and analysis of the self-regulation model. 

For dynamical equilibrium to be satisfied, the total 
pressure at the midplane must balance the gravitational 
weight of the overlying diffuse-ISM gas, P tot = W = 
(l/2)E(jigf(g z ). In different regimes, this pressure may 
be dominated by different terms (thermal, turbulent, or 
radiation), but each pressure term individually responds 
to the star formation rate. Where there is a substantial 
amount of atomic gas heated by stellar UV, balance of 
heating and cooling leads to an equilibrium thermal pres- 
sure Pth oc J uv oc Esfr (OML10, KKOll). Similarly, 
balancing turbulent driving associated with expanding 
radiative SN remnants (or other massive-star momen- 
tum sources) with dissipation on a vertical crossing time 
leads to an equilibrium turbulent pressure Pturb oc Esfr 
(Paper I, KKOll). In extremely high E regions, trapped 
reprocessed starlight provides a radiation pressure P ra d oc 
EEsf r that begins to comp ete with the turbulent pres- 
sure ([Thompson et alll20 05. Paper I). Putting these indi- 
vidual terms together, P tot = P th + Pturb + P-ad oc E S fr- 
Thus, under self-regulation the combined constraints of 
thermal, turbulent, radiative, and dynamical equilibrium 
imply that the star formation rate will naturally evolve 
to a level imposed by the vertical gravitational field, 
Esfr oc W. 

In mid- and outer-disk regions (generally where E ,$ 
100 Mq pc" 2 ), the warm (T ~ 10 4 K) ISM is space- 
filling and GMCs appear to be self-gravitating struc- 
tures that do not participate in the general vertical equi- 
librium. For this regime, OML10 show that the ther- 
mal/dynamical equilibrium theory is in good agreement 
with observations, with Esfr depending on both E and 
the stellar density p* of the disk (see also KKOll). For 
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outer disks, diffuse^ HI dominates and Esfr oc ^\fp* 
because the weight of the diffuse ISM is W oc S-^/pT 
in this regime. For mid-disks, gas is concentrated in 
gravitationally-bound clouds (observed as GMCs) which 
have relatively uniform column density, star formation 
efficiency, and other properties (probably as a result of 
internal feedback), such that Esfr oc E. 

In very dense regions where E ^ 100 M pc~ 2 , such 
as the Galactic center and ULIRGs, UV heating is not 
expected to play a strong role, and molecular gas is per- 
vasive rather than concentrated in effectively isolated 
GMCs. The transition to the "diffuse molecular" star- 
burst regime occurs where the pressure of the ISM as a 
whole exceeds the pressure of isolated, bound GMCs as 
found in the outer disk. For bound or virialized GMCs 
with surface density Eqmc = M/(ttR 2 ) that have a 
gravitational-to-kinetic energy ratio of 1 to 2, the inter- 
nal pressure is (0.5 — 1)GEq MC . The ISM as a whole 
must have midplane pressure (7r/2)GE 2 if equilibrium 
holds and gas dominates the gravity (see Equation [TO] 
below); from Paper I, this pressure is primarily tur- 
bulent, driven by star formation feedback. The tran- 
sition to the regime where molecular clouds lose their 
identity (and may be destroyed by externally-driven tur- 
bulence rather than internal feedback) therefore occurs 
when E > (0.5 - 0.8)E G mc ~ 50 - 100 M pc" 2 . 

In the starburst regime, the theory of Paper I suggests 
that SN play a key role in controlling the overall star 
formation rates because they dominate the momentum 
injection rate to the ISM0 Paper I presented the analyt- 
ical theory, compare d the star formation rates to obser- 
vations compiled in iGenzel et al.l (|2010f ) , and provided 
initial results from numerical models of SN-driven tur- 
bulent feedback in a cold ISM. According to the theory 
of Paper I, Esfr oc E 2 is expected for most starbursts 
(see Equation [TBI below) . Here, we extend Paper I to test 
the predictions from self-regulation over a wide range of 
galaxy and feedback parameters, using time-dependent 
numerical simulations. 

1.3. Simulations of Self- Regulation Due to Feedback in 
Starbursting Environments 

In this work, we model the evolution of a molecular 
dominated ISM using radial- vertical simulations, includ- 
ing a treatment for gas motion in the azimuthal direction. 
Using a large suite of hydrodynamic models, we focus on 
the role of SN driven feedback in the starburst regime, in- 
cluding its relationship to other disk characteristics such 
as the overall star formation rate, disk thickness, and 
midplane density. A key feature of these simulations is 
that the vertical dimension is well resolved, which is im- 
portant for accurately capturing the effect of turbulence 
on dis k thickness, as pointed out by iShettv fc Ostriken 
(2008). We test the sensitivity of the results to the as- 

3 We use the term "diffuse" to refer to spatially dispersed gas 
(both warm intercloud medium and cold cloudlets) that does not 
occur in gravitationally bound molecular clouds; see Section 2.2 of 
OML10. 

4 While thermal gas pressure from H II regions and radiation 
pressure are likely most important in destroying individual outer- 
galaxy GMCs containing embedded clusters (because of the time 
delay before supernovae), simple estimates suggest that for the ISM 
as a whole, the momentum input /stellar mass formed is dominated 
by supernovae (see Paper I). 



sumed input parameters, such as the efficiency of star for- 
mation in dense gas, and the momentum injected per unit 
stellar mass. Our analysis aims to understand the role 
of feedback-induced turbulence on the self-regulation of 
star formation in high (surface) density regions, where E 
> 100 Mq pc~ 2 , representative of the ISM in (U)LIRGs 
and galactic centers. 

This paper is organized as follows. The next section 
describes the relevant equations and our numerical meth- 
ods. Section 3 presents our model results, as well as a 
comparison between the simulations and the predictions 
from self-regulation theory. After a discussion we sum- 
marize our work in Section 4. 

2. NUMERICAL METHODS 

2.1. Basic Equations and Local Disk Model 

To model the evolution of the ISM in dense molecular 
disks, we solve the time-dependent hydrodynamic equa- 
tions, including self-gravity. The relevant equations are: 

g + V-(pv) = (1) 

|J+VVV = --VP-2«X V- V$ s +gext (2) 



dt 



V 2 $ 9 =47rGp 



(3) 



where p, v, P, and f2 are the volume density, velocity, 
pressure, and angular velocity of the rotating frame, re- 
spectively, and G is the gravitational constant. Since 
gas cools efficiently in high density molecular regions, we 
employ an isothermal equation of state, with constant 
sound speed c s — (P/p) 1 ^ 2 . We implement a static stel- 
lar gravitational field g ex t = — Q 2 zz (assumed to arise 
from a spherical bulge), which helps to concentrate gas 
to the midplane. The bulge potential is also responsible 
for the overall rotation of the gas with angular velocity 
Q. The time- varying self- gravitational potential due to 
the gas is $ g . 

The domain of our simulations is two-dimensional 
(2D), consisting of a radial and vertical (R, z) cross- 
section of a galactic disk, with extents Lr. and L z , respec- 
tively. Though the 2D simulations only treat x = R — Rq 
and z as independent variables, velocities in all three 
directions (including </>) are included. We also include 
Coriolis forces, with fi = flz constant (i.e. solid body 
rotation, for a constant-density bulge). We do not con- 
sider shear, as our focus is on galactic central regions, 
where the rotation curve is still rising. When dft / dR = 0, 
the tidal potential term in the rotating frame is zero 
and does not enter the momentum equation (this tidal 
term is nonzero in outer-disk regions where rotation is 
strongly sheared - see the right-hand side of Equation 
15ofKK011). 

Additionally, in our calculation of star formation rates, 
we implicitly consider the extent in the azimuthal di- 
rection (see Section I3.3|) . To model a local patch 
of the disk cross-section, we adopt periodic boundary 
conditions in R. In order to maintain a constant value 
of E throughout the simulation, we also adopt periodic 
boundary conditions in z. As we describe in Section |3~21 
we ensure that L z is sufficiently large in order to follow 
the complete evolution of the supernova shells, such that 



4 



Shetty & Ostriker 



the ISM scale height and star formation rate are con- 
verged. Simulating 2D (R, z) slices allows us to perform 
calculations with very high (sub-pc) spatial resolution, as 
well as to explore a wide range in parameter space (which 
may be used as a basis for future three dimensional [3D] 
simulations; initial tests show that similar results hold 
for 3D models). 

We numerically integrate the hydr odynamic Equation s 
H|) - © using the Athena code (jStone et all 120081) . 
Athena solves the partial differential equations using 
a single-step, directi onally unsplit Godunoy method in 
multiple dimensions (|Stone fc Gardinerll2009f ). We adopt 
piecewise-linear reconstruction and the HLLC Riemann 
solver. To solve the time-varying self-gravitational po- 
tential <£> ff , we employ a Fourier transform method 
with vacuum vertical boundary conditions and peri- 
odic horizontal boundary conditions, as described in 
iKoyama fc Ostrike r (2009b). We explore a range in Lr 
and L z , as well as the number of zones N-r and N z , in 
order to ensure that the results are not sensitive to the 
domain extent and that the features are well resolved 
numerically, as we discuss in Section [3.21 

2.2. Feedback Prescription and Model Parameters 

Equations ([IJ - ((3]) only describe the basic hydrody- 
namics, rotation, gas self-gravity, and the vertical po- 
tential. Our simulations also include an idealized model 
of momentum feedback produced by supernovae, which 
drives turbulence and disperses dense regions. This feed- 
back mechanism increases the total pressure, and limits 
collapse of the gaseous disk to only a small fraction of 
the densest material. 

Our method to identify regions that could form stars, 
and to apply momentum feedback tha t these stars would 
produ ce, is similar to that described in lShettv fc Ostriker] 
(2008). Here, we provide an overview of thi s algor ithm, 
and refer the reader to iShettv fc Ostrikeil (|2008t ) and 
KKOll for a more detailed description. 

We employ a statistical approach to determine host 
locations for the feedback events, and how much star 
formation is tallied (we do not remove gas from the do- 
main). Star formation can occur in a fraction of the re- 
gions where the number densities are greater than some 
chosen threshold density n t h- Thus, at every time-step 
each grid zone with n > n t h is identified. Next, the num- 
ber of massive stars (that can produce feedback) in zones 
with n > n t h is determined through a probability de- 
fined by two other user-chosen parameters, the "free-fall 
efficiency" (conversion of gas mass to stars per free-fall 
time), eg (nth), and the total mass in all stars formed per 
high mass star, to*. We then apply feedback instanta- 
neously, centered on those zones where high mass stars 
are determined to form (i.e. we omit time delays and 
spatial offsets in feedback, which more realistic models 
would take into account). The probability of a feedback 
event centered on a zone with n > nth in a given time- 
step At is thus 

p = Ate ff (rc t h)M c i ^ 
tff (n ai ) m t 

where M c \ is the mass of gas contained in the dense 
cloud in which the event originates, and the free fall time 
is tff(nth)= [37r/(32G/iTOj,7T, t h)] 1 ^ 2 ; here (i is the mean 



molecular weight and m p is the proton mass. For each 
massive star formed in a given time step, the total mass 
in stars formed is augmented by to* (Equation 21 of 
KKOll). 

After a zone is determined to host a supernova, a cir- 
cular region with chosen radius R s h is delineated. The 
density inside this region is reset to a uniform value (con- 
serving total mass), and velocities pointing away from 
the center are set such that the mean (spherical) mo- 
mentum injected per event is equal to p* (see Equation 
23ofKK011). 

In summary, there are five user-defined parameters re- 
quired to identify and implement feedback: n t h, eg (nth), 
Rsh, p*, and to*. We adopt to* = 100 M Q , which is de- 
rived from a iKroupal (|200lD IMF assuming supernovae 
result from stars with mass > 8 M Q . The chosen value 
of R s h also sets the effective azimuthal thickness L^— 
2R s h, which is used in setting M c \. The remaining three 
parameters, along with S, fl, c s , Lr, L Zl and the res- 
olution NrxN z complete the set of inputs for each nu- 
merical simulation. Table [1] lists the symbols and the 
corresponding description of the relevant model parame- 
ters and measured quantities we refer to throughout this 
paper. 

Our initial vertical density profile decreases as a Gaus- 
sian away from the midplane, such that the surface den- 
sity is S. We also include a sinusoidal perturbation along 
R, to seed gravitational instability. We have verified that 
our particular choice of initial conditions does not affect 
the later evolution in any way. As we demonstrate in the 
next section, by approximately one orbital time i or b= 
27r/0, the dynamic disk settles into a statistical steady 
state, such that the downward motions due to the verti- 
cal potential are countered by the upward motions due 
to feedback occurring near the midplane. 

2.3. Missing Physics 

The hydrodynamic models we consider are highly ide- 
alized, while in the real ISM a number of additional phys- 
ical processes may play a role. Cosmic rays, magnetic 
fields, and thermal radiation can contribute pressure, and 
can in principle affect self-regulation of star formation. 
The first two are, however, likely to be less important 
than the turbulent pressure if cosmic ray and magnetic 
scale heights are large compared to that of the neutral 
disk, and the last is likely important only if the gas sur- 
face density is extremely high (see Paper I). The ana- 
lytical model for self-regulated star formation in Paper 
I allows for feedback processes in addition to the turbu- 
lent driving considered here, and it will be interesting to 
explore these effects quantitatively in future simulations. 

As our simulations represent radial-vertical slices 
rather than full three-dimensional regions, we cannot 
study the detailed morphological structure of the ISM, 
such as filaments and the shapes of dense clouds. Three- 
dimensional simulations would be necessary to charac- 
terize the masses of clouds, and to make comparisons 
to structures as identified in po sition-position- y e locity 
molecular-line data cubes (e.g. iPichardo et all 120001: 
[ Ostriker et all 120011 : iGammie et alJ 120031 : IShettv et al.l 
|201Q[ ). Because the interior of vertically-expanding shells 
can be "filled" by gas moving horizontally from other 
azimuthal locations, the morphology in our present sim- 
ulations appears more "open" than it would in a fully 
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TABLE 1 
Symbols Employed 



Symbol 



Definition 



Simulation Parameters 



e fi{ n th) free-fall efficiency at the threshold density 

Ljt physical extent in radial dimension 

L z physical extent in vertical dimension 

m* total mass of stars per feedback event 

jVjt number of zones in radial dimension 

N z number of zones in vertical dimension 

nth threshold number density for feedback to occur 

Q angular velocity 

p» injected momentum per feedback event 

Fish radius of feedback event 

£ gas surface density 

t or t, orbital time 



Measured Quantities 




eff(rco) 


free-fall efficiency at midplane density 


fp 


turbulent dissipation parameter 


H 


gas disk thickness 


no 


gas number density at midplane 


^"drivc 


vertical momentum injection rate per unit area 


f-turb 


midplane turbulent pressure 


SsFR 


star formation rate 


tj„ 


vertical velocity dispersion 


Vz 


vertical velocity 


w 


vertical weight of the gas layer 


X 


contribution of bulge to vertical gravity, relative to gas self-gravity 



three dimensional model. We note, however, that three- 
dimensional simulations of self-regulated star formation 
in outer disks analogous to the radial-vertical models 
of KKOll give star formation rates that are quite con- 
sistent with those obtained using radial-vertical simula- 
tions. 

Because the primary focus of this work is on star for- 
mation in the molecule-dominated turbulent ISM, we 
have adopted the same (highly idealized) assumption 
of an isothermal medium that has been so fruitful in 
many of the first numeric al studies of turbulent mol ecu- 
lar clouds (see reviews b y IMac Low fc Kle sscn 20Q3 and 
iMcKee fc Ostrikerl [20071 ). In reality, the ISM has much 
more complex thermal and chemical structure, and a 
number of recent numerical studies have taken these 
into consideration. In particular, three-dimensional sim- 
ulations including detailed heating and cooling for ISM 
models with thermal super nova energy injection have 
recent l y been conducted b y Ide Avillez fc BreitschwerdtJ 
(|M^ : tlorm"g et all ([200^ ; ]HiiTet al.l (|2012T) . among oth- 
ers. Although most simulations including a hot ISM have 
focused on c o nditio ns similar to the Solar neighborhood, 
IJoung et al.l (|2009[ ) included a case with very high su- 
pernova rate, as would be expected for star formation 
rate ~ 1 M© kpc" 2 yr _1 , similar to the starburst regime 
we consider here. These recent multiphase simulations 
have not included self-gravity, however, and thus the su- 
pernovae rate is imposed as an input parameter rather 
than being modulated by the mass of gravitationally- 
collapsing gas. It will be quite interesting to include 
self-gravity and a feedback implementation together with 
multiphase heating and cooling to model self-regulated 
star formation more realistically. In particular, by com- 
parison with simulations that model supernovae by in- 



jecting thermal energy, it will be possible to assess the 
simple momentum injection model we adopt here to rep- 
resent turbulent driving in the neutral ISM by radiative 
supernova remnants. 

3. RESULTS 
3.1. Overview of Simulations 

We have explored a large range in simulation param- 
eters in order to develop a robust understanding of the 
effects of momentum feedback in high density, rotating 
disks. Table [5] lists the main simulations we consider 
here. We classify the simulations into five groups, based 
on the parameters which are varied. Column (1) indi- 
cates the name of each simulation, as well as the group 
to which it belongs. Columns (2) - (6) list the input 
values of surface density S, star formation efficiency per 
free-fall time at the threshold density eg (nth), momen- 
tum injected per supernova p», rotational speed fl, and 
orbital time i or b, respectively. The last two columns list 
the R and z dimensions of the simulation domain. Notice 
that some models are repeated in different Series: S100 
= PA3, S200 = E0.005 = 02, and E0.01 = PB3. We fur- 
ther note that although we have executed and analyzed 
well over 100 additional simulations, those listed in Ta- 
ble [5] span a sufficiently broad range of the parameters 
to highlight the major findings of our research. 

We have also explored variations in the other param- 
eters required to execute the simulations: m*, c s , nth, 
Rsh, £r, L z . As we discuss, m* always occurs as a ratio 
with in the relevant equations, so any variation in 
is equivalent to a corresponding variation in p*/m*. We 
thus fix m* to 100 Mq, and vary p*. We vary about 
the value expected for a supernova that has reached the 
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TABLE 2 

Input Parameters of Hydrodynamic Models 81 



(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


(8) 


Model 


S 


£ff(™th) 


P* 


n 




Lr 


L z 




(Mq pc- 2 ) 




(M Q km s" 1 ) 


(Myr" 1 ) 


(Myr) 


(PC) 


(PC) 


Series S 


(variation in S) 














S100 


100 


0.005 


3 x 10 s 


0.1 


62.8 


120 


240 


S200 


200 


0.005 


3xl0 5 


0.2 


31.4 


60 


120 


S400 


400 


0.005 


3xl0 5 


0.1 


15.7 


30 


60 


S800 


800 


0.005 


3xl0 5 


0.8 


7.9 


30 


60 


Series E 


(variation in CfF^th)) 














E0.005 


200 


0.005 


3 x 10 s 


0.2 


31.4 


60 


120 


E0.01 


200 


0.01 


3 x 10 s 


0.2 


31.4 


60 


120 


E0.025 


200 


0.025 


3 x 10 s 


0.2 


31.4 


60 


120 


E0.05 


200 


0.05 


3 x 10 s 


0.2 


31.4 


120 


240 


Series PA 


(variation in p*) 














PA1.5 


100 


0.005 


1.5 x 10 s 


0.1 


62.8 


120 


240 


PA3 


100 


0.005 


3 x 10 s 


0.1 


62.8 


120 


240 


PA6 


100 


0.005 


6 x 10 s 


0.1 


62.8 


120 


240 


PA9 


100 


0.005 


9 x 10 s 


0.1 


62.8 


120 


240 


Series PB 


(variation in p*) 














PB1.5 


200 


0.01 


1.5 x 10 s 


0.2 


31.4 


60 


120 


PB3 


200 


0.01 


3 x 10 s 


0.2 


31.4 


60 


120 


PB6 


200 


0.01 


6 x 10 5 


0.2 


31.4 


120 


240 


PB9 


200 


0.01 


9 x 10 s 


0.2 


31.4 


120 


240 


Series O 


(variation in Q) 














Ol 


200 


0.005 


3 x 10 s 


0.1 


62.8 


60 


120 


02 


200 


0.005 


3 x 10 s 


0.2 


31.4 


60 


120 


04 


200 


0.005 


3 x 10 s 


0.1 


15.7 


60 


120 


08 


200 


0.005 


3 x 10 s 


0.8 


7.9 


60 


120 



a All listed models have N R x N z = 512 X 1024 

shell formation stage (e.g. iBlondin et aTi ri998): 

\10 oi erg/ Vcm 6 / 

this is insensitive to the ambient density no and approx- 
imately linear in the supernova energy. In all the simu- 
lations, we set c s —2 km s . This value is larger than 
the sound speed of cold (T £ 100 K) gas. Such values 
are necessary because without magnetic fields, shocked 
gas would result in unrealistically high density regions, 
and thus lead to very small time-steps in the numerical 
simulations. Since turbulent motions still dominate, and 
to partly account for these (unmodeled) magnetic effects, 
we set c s =2 km s _1 . We have verified that provided c s 
is small compared to the turbulent velocity, the precise 
value does not significantly affect the results. We also 
note that for the analogous simulations of KKOll, ini- 
tial tests show that inclusion of magnetic fields do not 
significantly alter the results. For the remaining input 
parameters n t h, -Rsh, and box size Lr, L z , we discuss 
effects on the disk evolution in the subsequent sections. 

3.2. Box Size and Resolution Tests 

Before presenting the simulation results, we verify that 
the choice of domain size and the numerical resolution do 
not affect the outcome. Since we employ periodic bound- 
ary conditions, the extent in z must be large enough such 
that gas flow across the z boundary is unimportant. Gas 
leaving the (top or bottom) z boundary returns through 



zones. 

the opposite boundary; if outflow velocities were large, 
there would be a corresponding spurious compression of 
gas toward the midplane by the returning inflow. By 
constructing sufficiently large vertical domains, we en- 
sure that there is little mass and momentum flux through 
the boundaries^! In addition to the size of the domain, 
the physical resolution must be sufficiently high to en- 
sure that any gaseous structures that form, such as the 
high density clouds, are well resolved so th at the Truelove 
criterion is satisfied ()Truelove et al.|[l997f) . 

Figure (pj,) shows how the steady-state Esfr (defined 
in next subsection) depends on box size for model PB6, 
for a given physical resolution. When the ratio of L z to 
the disk thickness H (also defined below) is small, the 
midplane density is artificially enhanced (as described 
above), triggering more cloud collapse and subsequent 
supernova explosions. The star formation rate decreases 
as L z /H increases, with fewer shells passing through the 
boundary. At large L z /H, Esfr converges to a limiting 
value. 

The momentum fluxes (pv^) through the top and bot- 
tom boundaries are ~43% of the momentum flux in the 
disk midplane for the simulation with L z =20 pc. In the 
model with L z =320 pc, the boundary momentum flux 
is ;$ 0.5% of the value at the midplane. Similarly, we 
measure large time-averaged vertical mass flows (/o|w z |) 

5 In reality, hot gas produced by supernovae and high-altitude 
material accelerated by radiation forces may escape as a wind; the 
current simulations focus on cold, dense gas and do not include 
these effects. 
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Fig. 1. — a) The effect of box size L^/H on EgFR- Points show the mean Sgp^ from model PB6, but with different extents in L z . b) 
^SFR from S200 models with different resolutions (7V Z = 2x7Vr). In both panels, the points correspond to mean values of Esfr in 10 Myr 
bins, beyond 50 Myr from the start of the simulation. 



at the vertical boundaries in models with insufficient ex- 
tents0 For the model with L z =20 pc, the ratio of (p\v z \) 
at the vertical boundaries to that in the midplane is ^ 
0.6. When the vertical extent is sufficiently large, such as 
the model with L z —320 pc, this ratio is ^ 0.01. Based 
on a large number of tests of different models, we have 
found that the ratio of vertical mass flow and momentum 
flux in the vertical boundary to the corresponding value 
in the disk midplane is negligible when L z /H is large. 
Correspondingly, we find that Ssfr converges provided 
L z /H <L 6, so that for all simulations we choose a do- 
main size such that L z /H > 6, for the measured value of 
H. 

To ensure that the simulation results are independent 
of numerical resolution, we have executed a number of 
simulations with the same box size, initial conditions, 
and feedback parameters, but with varying and N z . 
Figure (Jib) shows the mean value of Ssfr for the fidu- 
cial model S200 with different resolutions, all with box 
size LrxL z = 60x120 pc 2 . Clearly, Ssfr converges 
to within 15% for all cases with dimension N^xN z > 
256x512. For N R xN z = 256x512, the physical reso- 
lution in model S200 is 0.23 pc; at our standard size 
NnxN z — 512x1024, the physical resolution is 0.12 pc. 
Our largest box is twice as large as that of model S200, 
with resolution 0.23 pc. At this resolution, the highest 
density at which the Truelove criterion (Aj/4 > L z /N z , 
for Jeans length Aj = c s [tt j '(Gp)] 1 ^ 2 ) is satisfied is ~ 10 5 
cm~ 3 , whereas typical cloud densities in our models are 
;$ 10 cm -3 . Thus, in order to explore a large range of 
parameters, and at the same time be confident that the 
simulations arc sufficiently well resolved, we will employ 
NrxN z = 512x1024 as the standard resolution. 

We have also explored the impact of the remaining 

6 The time-averaged true mass flux (pv z ) is zero at all heights. 



user defined parameters, R s h, n t h, and Lr. For ambient 
density of ~ 100—1000 cm~ 3 (similar to mean densities in 
our models) , supernova remnan ts become rad iative when 
their radii are a few pc (e.g. iDraind 1201 lL Equation 
39.21). We adopt a standard value of i? s h= 5 pc, and 
find similar simulation behavior for any other i? sn within 
a factor 2 of this value. We find that when n t h ~ 5000 
cm~ 3 , the evolution is not strongly dependent on the 
choice of n t h- 

Because we have periodic boundary conditions in the 
radial direction, the value of Lr does not affect the evo- 
lution provided that Lr ^ H and that the time-averaged 
gas distribution remains uniform and "disk-like" in the 
radial direction. However, for some conditions the value 
of the Toomre Q parameter (using the turbulent velocity 
dispersion) will be small enough that the combination of 
turbulence and rotational support is insufficient to pre- 
vent radial collapse under self-gravity for large radial do- 
mains (see Equation 29 of Paper I) . As discussed in Ap- 
pendix B of Paper I, the massive structures that form as 
a consequence of this collapse in real galaxies may poten- 
tially be dispersed by radiation pressure. However, in the 
current work we have not implemented radiation forces, 
so we consider only models that do not lead to overall 
radial collapse. For our fiducial model S200, the value of 
Q is less than unity, so that collapse ensues if we use a 
large radial domain. However, we have confirmed that 
if we increase f2 by a factor 1.5 (with other parameters 
as in the S200 model) such that Q > 1, a model with 
Ln = 120 pc is in all respects quite similar to the same 
model run with — 60pc; e.g. Ssfr differs by only 
~ 10%. Additionally, we considered a model similar to 
S200 but with fi=0.4 Myr" 1 andp*=6xl0 5 M Q km s" 1 
such that the Toomre parameter is in the stable regime. 
We find that the models with Lr between 50 and 100 pc 
achieve convergence in Ssfr to within ~ 5% 
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3.3. Model Evolution and Statistical Properties 

Figure [5] shows the volume densities of the four Series S 
models at two orbital times (2t rb) from the start of the 
simulation. As we discuss below, each model approaches 
a statistical equilibrium well before t m b'- the star forma- 
tion rate, vertical velocity dispersion, disk thickness, and 
other dynamical-state parameters all approach quasi- 
steady values. Numerous evolved SN shells are evident 
in Figure [2] One clear trend in Series S is that in mod- 
els with higher gas surface density, the gas is also more 
concentrated towards the midplane (note that the panels 
have different dimensions). 

The SN feedback events occur in the dense gas near the 
midplane, and are responsible for pushing gas to higher 
altitudes, as well as driving turbulence (both horizontal 
and vertical motions) and creating the filamentary fea- 
tures easily identifiable in Figure [2] Figure [3] shows the 
density of model S100 at t = 2.7t orh = 170 Myr, along 
with the velocities in the R, z plane. The large scale ve- 
locities are generally directed towards the midplane at 
this particular instant (althou gh at other times the over- 
all flow is expanding, e.g. see IWalters fc C^lifflOl . 

A close-up of two regions shows the detailed density 
and velocity structure. One region focuses on a patch in 
the midplane where a SN has recently exploded. The vec- 
tor field illustrates how gas within the SN shell is rapidly 
expanding away from the center of the bubble, even while 
surrounding gas is converging. The other close-up shows 
a region away from the midplane. The dense regions 
and filamentary structures evident here were created by 
interactions of gas driven by numerous earlier feedback 
events. Gas velocities near these dense structures devi- 
ate from the large scale converging flow towards the disk 
midplane. Feedback events thus influence gas motions 
far from their origin, driving turbulence throughout the 
simulation domain. 

In each model, the star formation rate Ssfr at time 
t is computed from the number of feedback events A^sn 
occurring over time interval Aibin centered on t. The 
contribution of mass to Ssfr is N^m*, where m* is the 
mass of all stars formed per star capable of undergoing a 
supernova. Since the star formation probability assumes 
an effective thickness of our simulation slice L^,= 2i? s h, 
the same effective thickness is used in computing the area 
of the domain projected on the horizontal plane, L^L^. 
Thus, Esfr over a given time interval is 



^sfr = 



TO*iVsN 

LnLcfjAtb 

h 



(6) 



Figure[4]shows the evolution of Ssfr, computed in bins 
of Atbm=20 Myr, as a function of time from all the Series 
S models. This value of Aibm is much larger than the 
vertical crossing time of each simulation, which is simply 
the thickness H of the disk divided by the characteristic 
vertical velocity v z , both of which are defined and an- 
alyzed below. We can thus be sure that the estimated 
SgFR is averaged over a sufficiently long time such that 
(on average) gas has cycled between the mid-disk z — 
and out-of-plane \z\ > locations numerous times. Fig- 
ure S] indicates that Ssfr saturates within 50 Myr, and 
as we discuss below in Section 13.41 the saturated value 
generally approaches the predictions from self-regulation. 

Two quantities describing the gas kinematics and disk 



structure are the velocity dispersion and disk thickness, 
respectively. To quantify the vertical motions, we com- 
pute the mass-weighted z-velocity dispersion through 



-,1/2 



(7) 



where the summation is taken over all zones in the sim- 
ulation. Figure ([5^) shows the evolution of the velocity 
dispersion for model S200. As does Ssfr, <?v also statis- 
tically converges, in this case to ~4.5 km s . Similarly, 
the mass- weighted disk thickness is defined as 



H = 



1/2 



(8) 



Higher surface (and volume) densities lead to thinner 
disks, as evident in Figure [21 Figure ((SJd) shows that 
H for model S200 saturates at ~ 9 pc. 

Given the vertical velocity dispersion and thickness, 
the vertical dynamical time is t vcr = H/a z ss 2 Myr for 
model S200. The measured quantities in Figure[5]are the 
mean values in 10 Myr bins, so that each bin corresponds 
to ~ 5 dynamical times. Again, this allows sufficient time 
for gas to cycle between the dense and diffuse phases. 

Another quantity of interest is eg (no), the efficiency 
of star formation per free-fall time, where iff (no) is eval- 
uated at the mean midplane density no. As discussed 
in Paper I, eg (no) represents the overall efficiency of star 
formation at the prevailing ISM conditions, and need not 
be the same as the value eg (nth) imposed to set the rate 
of star formation in very high density gas (see Sections 
3UandSl). 

Since the star formation rate can be directly measured 
through Equation ([6]), and iff (no) can be calculated from 
the (horizontally- and time- averaged) midplane density 
no measured in the simulations, the mean measured star 
formation efficiency is given by: 



e ff (n ) = £ SFR iff(no)/£. 



0) 



Figure ((5t) and ([5ji) respectively show the evolution of 
the midplane density, no, and the mean efficiency, eg (no), 
for model S200. As with a v and H, these quantities 
also saturate, with steady-state values no=385 cm~ 3 and 
e ff (n )=0.0041. 

Figures [6] and [7] show the mass- weighted density and 
velocity probability distribution functions (PDFs) for Sc- 
ries S models. The distributions show the average PDFs 
from times 50 Myr < t < 150 Myr, assessed in 5 Myr 
intervals. The simulations with higher surface densities 
produce PDFs which are systematically shifted towards 
larger volume densities. For model S800, the magni- 
tude of the (self-gravitational and external) potential 
strongly confines gas to the disk midplane, such that 
the disk thickness becomes comparable to our chosen 
value of the SN shell radius «5 pc (see Fig. [5]). The 
feedback events produce thin shells of shocked gas that 
have very high densities, which result in the high-density 
secondary peak. Yet, most of the mid-disk has den- 



sity 



5000 cm , corresponding to the main peak 



in Figure [BJi. Apart from the S800 model, these den- 
sity PDFs are all well represented as log-normals, as 
expected for highly compressible turbu l ent flows (e.g. 
Vazqucz-Scmadcni 1994; Klcssen| 120001 : lOstriker et al.l 
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Fig. 2. — Density of Series S models at times </i or b ~2. 

2001). The mass-weighted velocity PDFs are approxi- 
mately normal, but have more pronounced tails at both 
high and low velocities. The velocity PDFs do not show 
any significant differences among the Series S models, 
indicating that turbulent velocities are not strongly de- 
pendent on S (or Ssfr), a point we return to in Section 
EOl 

Table [3] summarizes the mean values Ssfr 5 & v H, uq 
and eff(no) for all models. Averages are based on bins 
of 20 Myr, starting at t = 50 Myr. The last two columns 
give f p and x> quantities relevant for the analytical ex- 
pressions derived in Paper I and discussed below. All 
models with entries listed for the measured parameters 
reach a steady state. However, Model 08 did not col- 
lapse to reach high densities; we believe this is because 
it is stabilized by a high rotation rate (see Sec. I3.4[) . In 
the following section, we compare the properties of each 
of the simulations presented in Table [3] to the analytical 
predictions from self-regulation derived in Paper I. 

3.4. Comparison with Predictions from Self- Regulation 

As discussed in Paper I, in self-regulated starburst re- 
gions, the vertical weight W of the molecular disk is ex- 
pected to be balanced mainly by turbulent pressure Pturb 
(unless the optical depth to IR exceeds ~ 16). Under this 
framework, turbulence is driven predominantly by feed- 
back from massive stars, so the momentum injection rate 
determines Pturb- The total upward momentum per unit 
time per unit area is then given by /pPdrivc for Pdrivc a 
fiducial momentum injection rate per unit area associ- 
ated with star formation and f p an order-unity dimen- 
sionless constant. Each of the gravitational, turbulent, 
and feedback momentum-injection fluxes may be mea- 



sured directly in the simulations through: 
ttGS 2 



W = 



-(i + x), 



D 2 
-fturb — P0 a vi 



1 P* 

4 m* 



J SFRi 



(10) 

(11) 
(12) 



respectively. As discussed below, \ accounts for the 
gravity of the stellar bulge relative to gas self-gravity, 
and is usually small. In choosing a fiducial value for 
we assume that radiative supernova shells dominate 
the momentum injection (see Equation [5J, but other 
terms could equally well be included in Equation (jl2[) . 
and we explore a range of p*. If the disk evolves to 
be turbulence-dominated and governed by star forma- 
tion self-regulation, then we should find that Pt U rb~ 

Pdrivc- W. 

Figure [S] shows the relationships of the measured mo- 
mentum fluxes Pturb and Pdrive with W. We compute 
Pturb in the simulations using midplane horizontal- and 
time-averages of (pv^). The turbulent and SN momen- 
tum fluxes are in excellent agreement with the vertical 
weight of the disk. For those models showing the largest 
deviation from the expectations in Figure ([5^), simula- 
tions PB1.5 and S800, there is only a factor of two dis- 
crepancy between P tur b and W. In Figure ([Hp), models 
S400 and S800 have the largest discrepancy between the 
predicted and measured momentum injection rates. For 
very strong gravity models, the disk thickness becomes 
comparable to the (imposed) radii of SN shells in our 
models. As a consequence, the disk can become "arti- 
ficially" thickened, because real feedback shells starting 
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Fig. 3. — Densities of model S100 at time t = 2.7t or b = 170 Myr. Vectors show radial-vertical velocities (y^ + v 2 ) 1 / 2 . White vectors 
displayed in the bottom left of each panel show the vector scale. Dots indicate locations where the velocity is < 3 km s — 1 . The large box 
is 120 pc X 240 pc, and the inset boxes are each 45 pc X 29 pc. 



at much smaller radii and conserving momentum might 
not expand as much. If shells reach larger sizes than 
their "natural" radii, the corresponding mean density 
and pressure would be somewhat lower than would be re- 
quired for self-regulated equilibrium. Overall, the general 
correspondence between Pturb, -Pdrivo, and W strongly 
supports the idea that the evolution of our ISM models 
reaches an equilibrium governed by star formation self- 
regulation. To further explore this premise, we now turn 
our attention to comparing other physical properties of 
the simulations with the predictions from self-regulation 



theory. 

We begin by providing an overview of the analyt- 
ical results expected under self-regulation. The star 
formation rate in equilibrium is obtained by equating 
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a The values listed in this Table are depicted in Figures l9ll4l where the la variations are 
also provided. 
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fturb with W (Equation 13 in Paper I): 

_ 2tt(1 + x) m*GY? 
fp P* 

= 0.092 M kpc" 2 yi' _1 



100 M© pc- 2 



/ p v 3000 kms " 



(13) 



The factor \ accounts for the gravitational potential due 
to the bulge (see Section 2 and 4 in Paper I), with 



X = 



2C 



1 + VI + 4C 



(14) 



Here, C « 0.66W 2 , where W = (J V VL/ '(ttGE) is a pa- 
rame ter analogous to the Toomre Q parameter ([Toomrd 
|1964|) . Using typical values, 
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Fig. 4. — The evolution of Esfr with time for Series S models 
with gas surface densities S= 100 (black circles), 200 (red dia- 
monds), 400 (blue triangles), and 800 (green squares) Mq pc -2 . 
Points show Ssfr m temporal bins of Atb; n =20 Myr, as computed 
from Equation ( f6l . 



VlOkms" 1 / 



n 



0.1 Myr" 



100 M Q pc" 



(15) 

thus C is typically small in our simulations. 

The parameter f p characterizes the magnitude of tur- 
bulent dissipation, with f p ^l for strong dissipation and 
f p ~ 2 for weak dissipation. The value of f p is defined 
by the ratio of Pturb and the fiducial vertical momentum 
flux injected by star formation, Pdrive (see Equations [Til 
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bins. 



fp — Pturb 

=w 



( p* 



\4m* 
P* 



J SFR 



4m* 



J SFR 



(16) 
(17) 



have verified that f p measured through Equation (|16l) 
provides similar values, since Pturb ~ W (Fig. [17]). Table 
[3] also shows that \ is measured to be rather small. 

By equating Pturb and W, the turbulent velocity dis- 
persion can be expressed as a relationship between the 
characteristic vertical acceleration under self-gravity, ~ 
a v /tff, and the mean gravitational field ~ t:GY<: 



where the second line assumes that dynamical equilib- 
rium also holds (see Fig. [5J as well as Fig. 8 of KKOll). 
Accordingly, in a self-regulated system, 



a v = ^=t s (n )Gi;(l + x) 1/2 



fp- 



= 0.92(1 + *) 

p*/m, 



J SFR 



O.1M kpc 



-2 yr 1 



lOOM pc" 



= 4.42 km s" 1 



100 cm- 



-1/2 



100 M pc" 



:(1 + X) 



1/2 



(19) 



3000km s- 



(18) 



this is a simply a re-arrangement of Equation ([IT 
We note that Pturb = Po^v is equivalent to 
(Tla v /2)(H/a v )~ 1 . With Scr„/2 the vertical momentum 
per unit area contained in each side of the disk, and H/a v 
the vertical crossing time, the relation Pturb ~ Pirivc, or 
f p ~ 1, thus implies that the disk's vertical momentum 
is replenished by feedback approximately once per dy- 
namical time. 

Using the value of SgFR measured using Equation ([6]) 
and x from Equation (|14[) . we can calculate f p in each 
simulation from Equation (|18[) . To obtain \ through 
Equations ([1^)) and (fT5)l . we use the measured value of 
a v . Table [3] provides the values of f p measured from the 
simulations in this way; all values are near unity. We 



Equation (|19p should hold for any disk-like system 
supported primarily by turbulence, independent of the 
source of that turbulence. 

The predicted velocity dispersion a v can also be ex- 
pressed in terms of eg (no), f P , and x as: 



<j v — 5.5 km s 



fp 



eff(n )\ / p*/m* 



(1 + X) 1/2 V °- 005 / \3000 kms- 1 

(20) 

(see Equation 22 of Paper I) ; Equation (|2U| follows from 
Equation (|19[) using the definitions of eft (no) and f p from 
Equations © and (TT51) . respectively. This form shows 
that if f p and eg (no) are approximately constant, then 
the velocity dispersion would be proportional to the mo- 
mentum/mass injected by star formation. 
Lastly, the predicted disk thickness H when dynamical 
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equilibrium holds is: 



1 



1 + x ttGS 
1 



= 74 pc- 



1 + X ^10 kms 



100 M Q pc- 



^l) 



(Note that the first equality in Equation 28 of Paper 
I contains a typo; the denominator should contain a E 
instead of a E 2 .) Using Equation (f2"0"|). this may be re- 
expressed as 



i? = 23pc 



fp ( eff(no) 



{i + xY 



0.005 



V3000 kms" 



lOOM pc- 



(22) 



When the definitions for eg (no) and f p (Equations [9] 
and [TBI are substituted into Equation (|22p . the result is 
E/(2p )- While Equation (|2"Tj) should hold independent 
of the source of turbulence, Equation (f22j) shows that H 
would scale inversely with E for self-regulated turbulent 
disks if f p and eg(no) remain approximately constant. 

Using the measured values of Esfr, Cu, H, and no in 
each simulation, we can test a number of aspects of the 
theory in Paper I. In particular, we can: (1) compare 
our measurements of Esfr to Equation (I13p to assess 
the combined (turbulent driving/dissipation and grav- 
ity/pressure) equilibrium and test whether f p ~ 1 is 
satisfied (for varying physical parameters E, p„, f2 and 
varying numerical parameter eg (nth)); (2) compare our 
measurements of a v to Equation (fit?)) to assess the bal- 
ance of turbulent pressure and weight, also comparing to 
Equation (|20|) to evaluate whether f p and eg (no) are ef- 
fectively constant (for varying parameters); (3) compare 
our measurements of H to Equation (]2 1 1) to assess dy- 
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predicted in Equation (fl~3l) (1.6 vs. 2), a larger adopted 
eff(^th) leads to a slightly steeper slope, so that our over- 
all results are generally consistent with Esfr °c E 2 (see 
Fig. 4 of Paper I). 

Equation (p~3]t indicates that Esfr under self- 
regulation is independent of the star formation efficiency 
in dense gas. Figure [9Jd indeed shows that the measured 
value of Esfr for Series E models is relatively insensitive 
to the chosen value of efr(nth)- Physically, this means 
that (within limits) the rate of star formation in dense 
gas does not affect the overall star formation rate av- 
eraged over large scales, because the amount of mass 
at high density simply adjusts until the feedback rate 
matches what is required to produce the needed turbu- 
lent pressure. 

From Equation (fT3|) . the star formation rate in self- 
regulated equilibrium should be inversely proportional 
to the input momentum per stellar mass p*/ m *: where 
p* is associated with high-mass stars and m* includes 



namical equilibrium, also comparing to Equation (|22[) to 
evaluate whether f p and eg (no) are effectively constant 
(for varying parameters). In addition, we can (4) use 
our measurements of Esfr and hq to compute a mea- 
sured eg (no ) via Equation (j9]) and explore whether there 
are any systematic dependencies on the physical or nu- 
merical parameters. 

Figure [9] shows the mean Esfr for all models after a 
steady state is reached (generally t > 50 Myr). The star 
formation rate is plotted against the main user-defined 
parameters varied between models from each series, a) 
E (Series S), b) eg (nth) (Series E), c) p* (Series PA and 
PB), and d) (Series O). The dashed lines in each panel 
show the predictions from self-regulation theory (Equa- 
tion ll3[) . for f p = 0.5 and 1.5, and with x=0. 

For Series S, Figured shows a remarkably good agree- 
ment between the measured star formation rate and the 
prediction for f p ~ 1. Although the increase of Esfr 
with E for Series S is slightly shallower than the power 
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all of the lower-mass stars proportionally (based on the 
IMF). Figure|9j: shows Ssfr as a function of for Series 
PA and PB models. An inverse proportionality between 
EgpR and p„ is evident, comparing favorably to the pre- 
diction from self-regulation. 

The rate of star formation in a self-regulated system 
is not expected to depend on the angular velocity f2, 
provided that angular momentum does not limit local 
collapse (i.e. on scales ^, H) and that the vertical stel- 
lar gravity is small compared to the vertical gas gravity 
(i.e. x ^ !)■ The independence of EgpR from is 
shown in Figure [9ji- We found, however, that when Q, is 
large enough - as in Model 08, angular momentum pre- 
vents clouds from collapsing to reach high densities (nth 
— 5000 cm~ 3 ), and we register it as non-star- forming. 
This model has Toomre wavelength At = 7r 2 GE/f2 2 « 
14 pc. This value of At is comparable to what would 
otherwise be the collapse scale, thereby stabilizing the 
ISM and preventing the formation of any clouds. 

Figure [10] shows how a v as measured in each simulation 
(Equation [7} compares to the expectation from vertical 
dynamical equilibrium (Equationll9j). There is generally 
a good correspondence between the predicted and mea- 
sured values. This comparison contains essentially the 
same information as in Figure (|BtO, and similar to the 
results there, the measured a v for a few models depart 
somewhat from the prediction. The greatest departure 
is for model S800, which is expected since the disk thick- 
ness approaches the numerically-imposed feedback shell 
size. 

Figure [TTJ shows a v measured in the simulations for 
each series. The dashed lines in each panel indicate the 
prediction from self-regulation given by Equation (|20[) . 
again with f p =0.5 and 1.5, along with x=0j an d using 
the mean value of eg (no) for each series. The predicted 
independence of a v from E and f2 is confirmed in Figures 
[TTk andlTTH. respectively. 



Figures [TTb-c indicate that the measured a v increases 
with eff(n t h) and p„, respectively. In Figure II lb . 
the lines correspond to Equation ([20]) with constant 
eff(no)=0.006. However, as discussed below, the mea- 
sured eff(no) increases with e ff (n t h) (also evident in Table 
[3]), implying from Equation (|2"0j) that o v should indeed 
increase with eff(n t h)- Similarly, Figure [TTh shows that 
the increase in <r v with p» is shallower than the linear 
relation indicated in Equation (|20l) for constant eg (no). 
This is due to the slight decrease in eg (no) withp,, which 
is further discussed below. 

Figure [12] compares the measured and predicted values 
of H, given by Equation ([5]) and (|2"TT) respectively. As 
with the velocity dispersion (Fig. [T0l) . the thickness is 
measured to be very similar to the predicted value. This 
agreement is indicative of vertical equilibrium between 
the weight due to gravity and turbulent pressure. 

Figure [T3l shows the measured thickness H from Equa- 
tion ©, compared to the prediction from Equation (|22p 
for constant f p and eg (no). Figures [T3h shows that H 
decreases with E less steeply than oc E _1 . Based on 
Equation (|22[) this is consistent with the systematic in- 
crease of f p with E for Series S (see Table [3]). Since the 
SN shell radius is chosen to be 5 pc in these numerical 
simulations, this places an effective lower limit on the 
disk thickness H 5 pc, which might be part of the 
reason for the shallow decrease of H with E. Similar to 
the case of a v , the increase of H with eg (nth) and in 
Figures [T3b-c can be fully accounted for by the variation 
of eg (no) with eg (nth) and p*, respectively. Figure [T3h 
shows that H is insensitive to f2, implying that neither 
rotation nor the external gravity of the bulge strongly 
affects the disk thickness, within the range shown. 

The measured value of eg (no), as computed through 
Equation ([9]) for each simulation, is shown in Figure 
1141 In general, there are only slight variations in eg (no) 
among the simulations; for Series S and O, eff(no) is ap- 
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(which also affects the vertical gravity). Error bars show the la deviation of the measured SgpR. 

proximately constant. As discussed above with regard 
to Ssfr, we interpret the weak dependence of eg (no) on 
eff(^th) as indicative of an adjustment in the mass of 
dense gas to meet the large-scale need for star formation 
feedback. This adjustment is possible because the dy- 
namical timescales decrease with increasing density and 
decreasing spatial scale. Other recent numerical stud- 
ies have also found that large-scale star formation rates 
are insensitive to user-defined parameters controlling star 
formation at small scales (see section l4~Tj) . Figure [T4"b 
demonstrates that eg (no) decreases somewhat with in- 
creasing . Potentially, this may be due to the increase 
of velocity dispersion with increasing p*, which renders 
a smaller fraction of gas eligible to collapse (see section 

ED) . 

In summary, based on our quantitative comparisons, 
the results from our numerical simulations show good 
agreement with the simple analytic theory of Paper I. 
Both vertical dynamical equilibrium and a balance be- 



tween turbulent driving and dissipation are satisfied. 
The dependence of Esfr, cr„, and H on the gas surface 
density £ and input momentum are similar to the pre- 
dicted behavior. In addition, the results are insensitive 
to the exact prescription for star formation in dense gas. 
The free parameter f p was introduced in Paper I to char- 
acterize the turbulent "yield" from momentum inputs by 
star formation, and our present simulations provide a nu- 
merical evaluation of f p . For our whole simulation suite, 
f p remains within ~ 50% of unity, the value for strong 
dissipation. 



4. DISCUSSION AND SUMMARY 

To investigate dynamics of the highly-turbulent, 
molecule-dominated ISM as found in (U)LIRGS and 
galactic centers, we have executed a suite of numerical 
simulations that incorporate feedback from star forma- 
tion. We demonstrate that in simulations reaching a 
steady state, many physical properties can be accounted 
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for by a simple theory of star formation self-regulation 
(Paper I). Namely, the turbulent pressure is driven by 
injected SN momentum, and dissipates within a vertical 
crossing time of the disk. The rate of star formation and 
momentum injection adjusts until the input rate of mo- 
mentum flux balances the vertical weight of the gaseous 
disk. 

4.1. Relationship to Previous Work 

Our numerical simulations of the ISM are similar 
in some respects to other recent modeling efforts that 
have included turbulent driving from localized feedback 
events, and our results are consistent with previous 
findings. In particular, we have found that the ve- 
locity dispersion a v is not strongly dependent on the 
the exact prescription for feedback as long as the mo- 
mentum (or energy) input is similar (|Dib et al.l 120061 : 



IShettv fc Ostrikerl2008HJoung et al.ll2009L KKOll). Ad- 
ditionally, we find that the overall star formation rate 
Esfr is not sensitive to the chosen value of eff(nt.h), in 
genera l ag reement with the con clusions of iDobbs et al.l 
(|2011f) and lHopkins et al.l (|2011| ) that the specific small- 
scale star formation prescription does not not strongly 
affect the resulting Esfr- Similar to previous ef- 
forts that have explored a large range of surface den- 
sities, our simulations here and in Paper I clearly 
demonstrat e a power law relationship between S and 
S SF R (e.g. iLi et al.l [20051 [200l iTasker k Bryan! [20061 



20081: iRobertson k Kravtsovl 120081: IShettv k Ostrikerl 
20081 IDobbs fe Pringlell2009t iKovama fc OstrikerT f2009a: 
Dobbs et al.ll2011b IHopkins et all 120111 ). Here, our nu- 



merical simulations are well resolved in the vertical di- 
rection, and we relate both the power law and coefficient 
of the Esfr. vs. E relationship to the requirements for 
equilibrium given in the self-regulation theory of Paper 
I. 

For the high-surface-density regime E 100 M© pc~ 2 



studied in this work, observations show that the Esfr 
vs. E relationship is steeper than linear. As discussed in 
Paper I, an accurate calibration of Xco , the ratio of gas 
mass to velocity-integrated CO intensity, is crucial for 
estimating S (and thus the exact power law of Esfr vs. 
E) from CO observations. Recent theoretical efforts have 
advanced our un derstanding of Xgo , in bo t h Milky- Way 
like G MCs (e.g. IGlover fc Mac Low! l20ll IShettv et all 
l2011allbT ) , as well as large scale galaxies and merger sys- 
tems (e.g. iNarayanan et al.ll2lTll"[ l2012f ). These models 
use a combination of numerical hydrodynamic simula- 
tions and radiation transfer to assess environmental de- 
pendencies of Xco- 

If gas dominates the vertical gravity, the theory of 
Paper I results in a power-law relationship Esfrcx E 2 
(Equation I13j) : the numerical simulations presented in 
Paper I and here (Fig. [9]) support this model. As 
demonstrated in Paper I, employing a continuously vary- 
ing Xco indeed shows Esfroc E 2 for a sample of 
ULIRGs and the Galac tic center (iGenzel et al.l 120101 : 
lYusef-Zadeh et al.l I2009D . INarayanan et al.l (|2012F in- 
vestigated the relationship between Xco and the ve- 
locity integrated CO (J = 1 — 0) brightness tem- 
perature Wqo m a large compilation of low- and 
high-z galaxies. Applying the model-based calibra- 
tion X co oc W^co' 3 ' INarayanan et al.l (|2012f ) found 
that E SFR ~ O.lM Q kpc- 2 yr- 1 (E/lOOM pc- 2 ) 1 - 95 , in 
agreement with the Paper I prediction (Equation ll3l here. 
with f p ~l and \ "C 1). 

The self-regulation theory of Paper I has a number of 
similarities to and differences from the star formation 
model in the high-surface-density molecule-dominated 
regime (E 100 Mq pc~ 2 ) proposed b y Krumholz and 
coworkers (|Krumholz fc McKed 120051 : iKrumholz et al.l 
120091 120121 ) . Both models rely on the role of supersonic 
turbulence. In Krumholz et al, the specific star forma- 
tion rate is characterized in terms of an efficiency per 
free-fall time at the mean density (essentially es(n,o)), 
where the mean density (w hich sets tff(no)) depends on 
E and the turbulence level. IKrumholz fc McKed (|2005[ ) 
argued that this efficiency should depend on the fraction 
of gas at pressures higher than the mean turbulent pres- 
sure, and pointed out that for log-normal density PDFs, 
this fraction depends only weakly on the Mach number 
(oc Ai~ 3 ) and is predicted to be ~ 0.01, consistent wit h 
observations of molecular gas ([Krumholz fc Taiil 120071) . 
Krumholz et al. do not, however, directly address the 
origin of turbulence - i.e. what sets a v . Rather, they 
adopt the assumption that Toomre Q (and therefore 
W) is order- unity, so that a v ~ TrGE/fi, and adopt 
an empirically-motivated relation Q cx E (so that 
a v cx E - 5 ) to obtain E SFR cx e ff E/i ff cx E 2 /^- 3 oc E 1 - 3 . 

Although the star formation rate in the current the- 
ory can also be characterized in terms of the velocity 
dispersion and es(no) (see Equation 21 of Paper I), the 
fundamental relationship is instead Equation () 13[) - This 
expression connects the star formation rate to the weight 
of the ISM (or equilibrium midplane pressure) and to the 
momentum/mass injected by feedback (p*/m»), yield- 
ing Esfr ~ 27rGE 2 (p*/m»)~ 1 . Equating this relation 

to E SFR = eff(no)E/tfif(no)« (4e ff (n a )/y r 3)GY, 2 /a v then 
leads to a proportionality between the velocity dispersion 
and both eg (no) and p»/m* (see Equation |20|) . Here, 
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we use numerical simulations to evaluate eg (no) (see 
Fig. Q3] and Table [3]), fin ding values in the range ~ 
0.00 5-0.01, consisten t withlKrumholz fc McICei (2005) 
and iKrumholz fc Tanl (|2007l ). We find, however, that 
the velocity dispersion is essentially independent of £ 
(see Fig. ITTa). which differs from the a v oc E 5 relation 
adopted by Krumholz et al. We note that due to lack 
of resolution, the molecular velocity dispersion on scales 
comparable to the disk thickness is difficult to obtain 
with observations, although this situation will improve 
with ALMA. 

4.2. Summary of Results 

We have conducted a suite of simulations in which we 
independently varied the gas surface density S, the input 
momentum per high mass star p*, the angular rotation 
rate of the gas Q, and the efficiency of star formation in 
very dense gas, eg("-th)- For each simulation, we mea- 
sured the star formation rate and ISM properties after a 
statistical steady state developed, and compared to the 



predictions of Paper I. Our main results are as follows: 

[1] For essentially all models, we find excellent corre- 
spondence between the turbulent momentum flux Pturb, 
the vertical weight of the gaseous disk W, and the vertical 
momentum injection rate per area Pdrivo associated with 
feedback (Fig. |5J|. The result that P turb ~ W w P d rive 
strongly supports the idea that the combined ISM/star 
formation system in starburst disks can be self-regulated, 
as described in Paper I. 

[2] Our results (Fig. show that Esfr is essen- 
tially independent of and eff(nth)) whereas Esfr in- 
creases for higher E and decreases for higher p* fol- 
lowing the expectations of self-regulation theory. The 
result that the large-scale Esfr is independent of the 
star formation rate in dense gas means that the pro- 
cesses with the longest timescales (associated with the 
largest spatial scales) are what controls the overall star 
formation rate. Physically this makes sense: gas that 
reaches high density collapses rapidly, but the (slower) 
rate at which this dense gas is resupplied by lower-density 
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Fig. 12. — The measured and predicted values of the di sk t hick- 
ness H for all models, as defined in Equations J5J and 12H . re- 
spectively. The agreement between the measured and predicted 
thickness indicates vertical equilibrium between the weight due to 
gravity and SN driven turbulent pressure. 

gas depends on larger-scale dynamics. As noted in Pa- 
per I (see also iNarayanan et al.l I2012T ). the prediction 
Ssfr ~ O.lMekpc^yr-^S/lOOMgpc- 2 ) 2 of Equa- 
tion ([T3|) also agrees with observations provided that 
Xqo decreases modestly with increasing £ (or Wqo)- 

[3] We find that the star formation efficiency per free- 
fall time at the mean midplane density, eff(^o)j is in- 
dependent of S, i7, and £fi(«-th), and decreases only 
slightly with increasing p*. The resulting erf(no)~ 
0.005 — 0.01 is similar to the theoretical estimates 
for tur bulent, self-gravitating gas at high Mach num- 
ber of iKrumholz fc McKed (|2005l ). while being some- 
what lower than the numerical estimates (from tur- 
bu lent simulations with perio dic boundary conditions) 
of iPadoan fc Nordlundl ([201 lh . Measured ratios of the 
stellar-to-gas content in nearb y molecular clouds are 
~ 0.03—0.06 (jEvans et a l. 2009), which would imply sim- 
ilar eg (rig) to our results if the cloud ages are several free- 
fall times. Gas at more extreme condit ions in ULIRGs 
is als o observed to have es(n a )^ 0.01 ([Krumholz fc Tanl 
12001 . 

[4] The vertical velocity dispersions in our models are 
in the range a v G 3 — 12 km s _1 for momentum per high 
mass star in the range € 1.5 — 9 x 10 5 M Q km s -1 . 
Similar to previous results, we find that a v is relatively 
independent of Esfr and also VI (see Fig. [TT]). The in- 
crease of v z with p* is shallower than linear, due to the 
decrease of eg (no) with increasing p» (see Equation [20]) . 



The agreement (Figs. [TU] and [T^]) between a v and H 
measured in the simulations and the respective predic- 
tions of Equations (fl9]) and (f2"Tj) shows that dynamical 
equilibrium between gravity and turbulent pressure is es- 
tablished. The disk thicknesses in our models are quite 
low (H ~ 4 — 40 pc) , increasing as increases and de- 
creasing as E decreases. 

[5] The densities and velocities in our models fol- 
low approximately log-normal and normal distributions, 
respectively (Figs. El [7|). These forms are a natu- 
ral consequence of supersonic isothermal turbulent flows 
(as discussed bv[Vazauez-Semadenll Fl994; K l essen! [20001 : 
lOstriker et~aT]|1999ll200lD . The log-normal density dis- 
tribution is a key feature invoked in various models 
of what sets the star formation efficiency in turbulent 
systems (IKrumholz fc McKe e 2005; Padoan fc Nordlundl 
mm iHennebelle fc Chabrierll2011h . 

A natural extension of the simulations presented here is 
to include the third dimension. High resolution 3D sim- 
ulations will allow for detailed morphological and kine- 
matic studies of the molecular ISM in starburst regions. 
Such simulations will also more accurately measure the 
parameter f p relating the turbulent pressure to the mo- 
mentum flux injected by feedback (Equation [16]). Fur- 
ther, more realistic modeling of the ISM should incor- 
porate a variety of feedback mechanisms and additional 
physics, including stellar winds and radiation, and heat- 
ing and cooling to follow cold, warm, and hot phases 
rather than an isothermal equation of state to follow just 
the cold gas. By combining with radiative transfer calcu- 
lations, such simulations will enable detailed comparison 
of feedback-regulated disks with observations of the ISM 
in starburst environments. 
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